Study of the
Sequence of Persistence Diagrams
Some simple examples
Vincent Runge
29/10/2023
library(tdaTS) #our package
library(plotly) #3D plotting
library(fields) #image.plot function1 data generators
The simple 3D plot function using function plot_ly
function from the plotly package.
myplot_3D <- function(data)
{
minD <- min(data$x, data$y)
maxD <- max(data$x, data$y)
M <- max(minD, maxD)
resScale <- rbind(data,c(0,minD,minD),c(0,minD,maxD),c(0, maxD,minD),c(0,maxD,maxD)) ### add 4 points in the corners for scaling x and y axes
plot_ly(resScale, x = resScale$x, y = resScale$y, z = resScale$t,
type="scatter3d", mode = "markers", color = resScale$t)
}We propose, as a starting point, the following four data generators. You give an example for each of the generators.
1.1 Segment to circle (isometric)
We close a segment into a circle, maintaining the length of the curve (2 pi)
data <- segment_to_circle(n = 1000,
change = 0.5,
time_sampling = "unif")
# an example with additional noise
data$x <- data$x + rnorm(n= 1000, mean = 0, sd = 0.2)
data$y <- data$y + rnorm(n = 1000, mean = 0, sd = 0.2)
######## plot the data
myplot_3D(data)1.2 Circle extinction
The circle is smashed into a segment equal to the diameter of the initial circle.
data <- circle_extinction(n = 1000,
change = 0.5,
time_sampling = "unif")
myplot_3D(data)1.3 Segment to Segments
we split a segment at the changepoint in two at the center point with
regular speed with final gap = 0.1 here.
data <- segment_to_segments(n = 1000,
change = 0.5,
gap = 0.1)
myplot_3D(data)1.4 Circle move and distortion
An example with no change. We deform the shape in a ellipse, both rotationg and shifting it.
data <- circle_move_distortion(n = 10000,
rotation = 2,
X_rate = 0.5,
Y_rate = 3)
myplot_3D(data)In most example, there is the possibility to use parameters
X_rate and/or Y_rate to extend or contract the
shape and then to change the density of points by level (not their
number)
data <- segment_to_circle(n = 1000,
change = 0.5,
X_rate = 0.5,
Y_rate = 3)
myplot_3D(data)2 Sequence of Persistence Diagrams
In following sections, we will filter the Diagrams using 3 possible thresholds:
birth <- 0.15
death <- 0.17
diag <- 0.1The birth threshold is a vertical threshold
The death threshold is a horizontal threshold
The diag threshold is a diagonal threshold
The objective would be to remove the noise part of each of the persistence diagrams.
We still need to find the right approach to filter the noise (considering filtering by ratio death/birth maybe?)
2.1 One example segment to Circle
In 3 steps
result with no noise
result with additional noise
segment comparison with matrix of cross Wasserstein distance using dimension 1 elements
set.seed(7)
data <- segment_to_circle(n = 1000,
change = 0.5,
time_sampling = "unif")
Plot_All_Persistence_Diagrams(data,
birth = birth,
death = death,
diagonal = diag,
nb_levels = 20)
data$x <- data$x + rnorm(n= 1000, mean = 0, sd = 0.1)
data$y <- data$y + rnorm(n = 1000, mean = 0, sd = 0.1)
Plot_All_Persistence_Diagrams(data,
birth = birth,
death = death,
diagonal = diag,
nb_levels = 20)
nb_levels <- 20
all_PD <- Generate_All_Persistence_Diagrams(data,
nb_levels = nb_levels)
################
#### thresholding
################
all_filter <- remove_noisy_points(all_PD,
birth = birth,
death = death,
diagonal = diag,
infinity = FALSE)
################
#### distances: Wasserstein
################
mydim <- 1 #choice for dimension in PD
t(all_filter[,1:2]) %>% as.matrix()## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## time 0 4 6 6 9 10 11 13 14 15 16 17 19
## dimension 0 0 0 0 1 1 1 1 1 1 1 1 1D_wasserstein_2by2 <- distances_Persistence_Diagrams(all_filter, type = "2by2", dimension = mydim)
par(mfrow = c(1,1))
D_wasserstein_2by2 <- t(D_wasserstein_2by2)
image.plot(D_wasserstein_2by2[,nrow(D_wasserstein_2by2):1], axes = F, col = grey(seq(0,1, length = 256)))
axis(1, at=seq(0,1,length.out = nb_levels), labels= 1:nb_levels)
axis(2, at=seq(0,1,length.out = nb_levels), labels= nb_levels:1)
segments(0,0.5,1,0.5, col = 2, lwd = 4)
segments(0.5,0,0.5,1, col = 2, lwd = 4)
2.2 Example circle_extinction
set.seed(7)
data <- circle_extinction(n = 2000,
change = 0.5,
time_sampling = "unif")
Plot_All_Persistence_Diagrams(data,
birth = birth,
death = death,
diagonal = diag,
nb_levels = 20)
2.3 Example segment_to_segments
data <- segment_to_segments(n = 1000,
change = 0.5,
gap = 0.5)
Plot_All_Persistence_Diagrams(data,
birth = birth,
death = death,
diagonal = diag,
nb_levels = 20)
2.4 Example circle_move_distortion
data <- circle_move_distortion(n = 10000,
rotation = 2,
X_rate = 5,
Y_rate = 1, time_sampling = "discrete",
nb_levels = 20)
data$x <- data$x + rnorm(n= 10000, mean = 0, sd = 0.1)
data$y <- data$y + rnorm(n = 10000, mean = 0, sd = 0.1)
Plot_All_Persistence_Diagrams(data,
birth = birth,
death = death,
diagonal = diag,
nb_levels = 20)
nb_levels <- 20
all_PD <- Generate_All_Persistence_Diagrams(data,
nb_levels = nb_levels)
################
#### thresholding
################
all_filter <- remove_noisy_points(all_PD,
birth = birth,
death = death,
diagonal = diag,
infinity = FALSE)
################
#### distances: Wasserstein
################
mydim <- 1 #choice for dimension in PD
t(all_filter[,1:2]) %>% as.matrix()## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## time 0 1 2 3 4 5 6 7 8 9 10 11 12
## dimension 1 1 1 1 1 1 1 1 1 1 1 1 1
## [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## time 13 14 15 16 17 18 19
## dimension 1 1 1 1 1 1 1D_wasserstein_2by2 <- distances_Persistence_Diagrams(all_filter, type = "2by2", dimension = mydim)
par(mfrow = c(1,1))
D_wasserstein_2by2 <- t(D_wasserstein_2by2)
image.plot(D_wasserstein_2by2[,nrow(D_wasserstein_2by2):1], axes = F, col = grey(seq(0,1, length = 256)))
axis(1, at=seq(0,1,length.out = nb_levels), labels= 1:nb_levels)
axis(2, at=seq(0,1,length.out = nb_levels), labels= nb_levels:1)
segments(0,0.5,1,0.5, col = 2, lwd = 4)
segments(0.5,0,0.5,1, col = 2, lwd = 4)
3 Removing the noise structure
Need to think…
Need to look at the litterature
Threshold diagonal ? death over birth ratio ?
3.1 Multiple circles in the 2d plane
library(TDA)
circle <- function(n = 1000, center = c(0,0), radius = 1, noise = 0.5)
{
pos <- runif(n, min = 0, max = 2*pi)
x <- center[1] + radius*cos(pos) + rnorm(n, mean = 0, sd = noise)
y <- center[2] + radius*sin(pos) + rnorm(n, mean = 0, sd = noise)
res <- matrix(c(x,y), nrow = n, byrow = FALSE)
return(res)
}
data <- circle()
DiagAlphaCmplx <- alphaComplexDiag(X = data ,
library = c("GUDHI", "Dionysus"),
location = TRUE,
printProgress = FALSE)
plot(DiagAlphaCmplx$diagram, col = 1 + cumsum(DiagAlphaCmplx$diagram[, 1]))
plot(data, col = 1,xaxt="n", yaxt="n",xlab="", ylab="", asp = 1)
one <- which(DiagAlphaCmplx[["diagram"]][, 1] == 1)
for (i in seq(along = one))
{
for (j in seq_len(dim(DiagAlphaCmplx[["cycleLocation"]][[one[i]]])[1]))
{
lines(DiagAlphaCmplx[["cycleLocation"]][[one[i]]][j, , ] + rnorm(n = 4,sd = 0.0), pch = 19, cex = 1, col = i + 1)
}
}